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We continue the study of collisionless systems governed by additive interparticle 
forces by focusing on the influence of the force exponent a on radial orbital anisotropy. In 
this preparatory work we construct the radially anisotropic Osipkov-Merritt phase-space 
distribution functions for self-consistent spherical Hernquist models with forces and 
1 ^ a < 3. The resulting systems are isotropic at the center and increasingly dominated 
by radial orbits at radii larger than the anisotropy radius Va- For radially anisotropic 
models we determine the minimum value of the anisotropy radius Vac as a function of a 
for phase-space consistency (such that the phase-space distribution function is nowhere 
negative for Va ^ Vac)- We find that r^c decreases for decreasing a, and that the amount 
of kinetic energy that can be stored in the radial direction relative to that stored in 
the tangential directions for marginally consistent models increases for decreasing a. In 
particular, we find that isotropic systems are consistent in the explored range of a. By 
means of direct IV-body simulations we finally verify that the isotropic systems are also 
stable. 


1. Introduction 


It is well known that spherically symmetric, self-gravitating collisionless equilibrium 
systems in which a signihcant fraction of the kinetic energy is stored in low angular 
momentum orbits are dynamically unstable. The associated instability is known as Radial 
Orbit Instability (hereafter ROI, see e.g. Fridman & Polyachenko 1984; Bertin 2014; 
Binney & Tremaine 2008; for a recent discussion see also Polyachenko & Shukhman 
2015 and references therein). Usually, the amount of radial anisotropy in a spherical 
system is quantified by introducing the Fridman-Polyachenko-Shukhman parameter 


2Kr 

where the radial and tangential kinetic energies are given respectively by 


(I.l) 


Kr = 27r 


p(r)(T^(r)r^dr, Kt = 27 r 


p{r)af{r)r‘^dr. 


( 1 . 2 ) 


In the expressions above p is the system density, and and ct/ are the radial and tan¬ 
gential phase-space averaged square velocity components (see Polyachenko & Shukhman 
1981; Polyachenko 1992): in particular, for isotropic systems ^ = 1. In the case of 


f Email address for correspondence: pierfrancesco.dicintio@unifi.it 
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Newtonian force numerical simulations show that the ROI typically occurs for ^ > 1.7, 
but this should be considered a fiducial value, as it is also well known that the critical 
value of ^ above which the given system is unstable depends on the specific phase-space 
structure of the equilibrium configuration under study (see, e.g. Merritt & Aguilar 1985; 
Berlin & Stiavelli 1989; Saha 1991; Berlin et al. 1994; Meza & Zamorano 1997; Nipoti, 
Londrillo & Ciotti 2002; Barnes, Lanzel & Williams 2009). 

Due to its relevant astrophysical consequences, and its intrinsic physical interest for the 
understanding of collisionless systems, the ROI has been extensively studied in Newtonian 
gravity (see e.g. Palmer & Papaloizou 1987; Gajda et al. 2015, and references therein). 
For example, Nipoti, Londrillo & Ciotti (2002) investigated the implications of the ROI 
for the Fundamental Plane of elliptical galaxies, excluding the possibility that the so- 
called tilt of the Fundamental Plane is just due to a systematic increase of radial orbits 
with galaxy mass (see also Ciotti 1997, Ciotti & Lanzoni 1997). In a different context 
(Ciotti & Pellegrini 2004) it has been shown that the constraints on the maximum 
amount of radial anisotropy that can be sustained by a stable stellar system can be 
used to dismiss some models of mass distribution in elliptical galaxies that have been 
obtained by using X-ray data under the assumption of hydrostatic equilibrium of the 
hot interstellar medium. The ROI also plays a role in the course of Violent Relaxation 
(Lynden-Bell 1967), i.e., the collapse and virialization of cold self-gravitating systems 
(e.g. see van Albada 1982; Londrillo, Messina & Stiavelli 1991; Nipoti, Londrillo & Ciotti 
2006), because such collapses are dominated in their last stages by significant radial 
motions, and so ROI contributes to establish the dynamical and structural properties 
of the quasi-relaxed final states (e.g. Merritt & Aguilar 1985; Trent! et al. 2005; Sylos- 
Labini et al. 2015). 

From a more general point of view, it is now known that the ROI is not a specific prop¬ 
erty of Newtonian gravity. For example Nipoti, Ciotti & Londrillo (2011) have studied the 
ROI in Modified Newtonian Dynamics (hereafter MOND; Milgrom 1983; Bekenstein & 
Milgrom 1984), comparing anisotropic MOND systems with their Equivalent Newtonian 
Systems with dark matter, (i.e., systems in which the phase-space properties and total 
gravitational potential are the same as in the corresponding MOND systems). Numerical 
simulations showed that MOND systems are more prone to develop ROI than their ENSs. 
At the same time, MOND systems are able to support a larger amount of kinetic energy 
stored in radial orbits than one-component Newtonian systems with the same density 
distribution. 

Some natural questions then arise. What is the dependence of the ROI on the specific 
force law? How much the long- or short-range nature of the interparticle force influences 
the ROI? Are there any aspects of the ROI that are peculiar of the Newtonian force 
(which stands out among power-law forces for several unique mathematical properties)? 
In fact, though several features of the dynamical evolution of collisionless systems have 
found to be quite independent of the specific force law, some other aspects show a clear 
dependence on it. Eor example, in the low-acceleration regime, collisionless relaxation 
is much less efficient in MOND than in Newtonian gravityf (Ciotti, Nipoti & Londrillo 
2007; Nipoti, Londrillo & Ciotti 2007a,b), and a similar result is found for power-law 
forces with a < 2 (Di Cintio & Ciotti 2011). 

In particular, Di Cintio, Ciotti & Nipoti (2013) (see also Di Cintio & Ciotti 2011) 
studied by means of direct V—body simulations the dissipationless collapses of cold and 
spherical systems of particles interacting via additive interparticle r““ forces. The main 

t Note however that in MOND collisional relaxation seems to be faster than in Newtonian 
gravity. See e.g. Ciotti & Binney (2004); Nipoti et al. (2008). 
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results of these works are that, almost independently of a, the considerably radially 
anisotropic final states of cold collapses have surface density profiles following the Sersic 
(1968) law, and differential energy distributions well described by an exponential over 
a large range of energies. At fixed initial virial ratio however, a non-monotonic trend of 
the Sersic index m (the parameter determining the shape of the density profile, see e.g. 
Ciotti 1991, Ciotti & Berlin 1999) with the force index is found. Remarkably, the a = 1 
case (corresponding qualitatively to the force in the deep-MOND regime), behaves very 
similarly to MOND, although this last theory is non-additive (the MOND field equation 
involves the p-Laplace operator). 

Prompted by these findings and by the above questions, here we set the stage for 
the study of ROl in the case of power-law forces, by constructing a family of radially 
anisotropic Hernquist (1990) models with additive r““ forces, and restricting to the 
range 1 ^ a < 3. Such interval contains the Newtonian a = 2 case, explores ’’short- 
range” forces (a > 2) and also long-range forces (a < 2), down to the MOND-like case 
a = 1. Smaller values of a are excluded by this preparatory investigation (we recall that 
the a = —1 case corresponds to the harmonic force, for which no relaxation takes place, 
see e.g. Lynden-Bell & Lynden-Bell 1982). Power-law forces with a ^ 3 are also excluded 
by the present analysis, due to their very peculiar properties (for example, self-gravitating 
systems with a = 3 can be virialized only for zero total energy). We notice that the study 
of the dynamics under the effect of power-law forces is by no means new: very important 
examples can be traced back to Newton’s Principia (e.g. see Chandrasekhar 1995). 

This paper is structured as follows. In Section 2, we construct isotropic and radially 
anisotropic self gravitating spherical Hernquist models with r““ forces, and we study 
their phase-space consistency and the properties of their velocity dispersion profiles. In 
Section 3 we presents numerical results on the stability of the isotropic systems. Section 
4 summarizes. 


2. Radially anisotropic equilibria 

The Hrst step for the realization of the A^—body simulations of ROI for systems ruled 
by generalized power-law forces is the construction of their equilibrium phase-space 
distribution function /, in a framework allowing for tunable amount of radial anisotropy 
(i.e. different values of f). It is well known that the amount of radial orbits that can be 
imposed to a physically acceptable equilibrium stellar system is limited by phase-space 
consistency (i.e. by the positivity of the associated /): stability can be studied only for 
consistent systems. Therefore, before embarking in the study of the stability of a system 
it is fundamental to be reassured that its equilibrium initial condition actually exists. For 
these reasons in the next Section we proceed to construct the phase-space distribution 
function for the widely used Hernquist (1990) profilef, under the effect of the r~°‘ force. 

2.1. The phase-space distribution function 

We begin by considering the problem of the construction of the phase-space distribution 
function / for a generic equilibrium spherical system of given density profile p(r), under 
the action of the force with 1 ^ a < 3. In principle, different velocity anisotropy 
profiles can be imposed on a given density proHle. In this work, in order to build 
systems with a tunable degree of radial anisotropy, we adopt widely used Osipkov-Merritt 

f A Referee pointed out to us that the Hernquist (1990) profile has been found earlier as 
limiting case of a generalized isochrone model (Kuzmin & Malasidze 1970; Kuzmin & Veltmann 
1973, see also Osipkov 1978) 
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extension of the classical Eddington (1916) inversion: 

1 d dpa dF 1 /-Q^up 

^ ~ y87r2 dQ Jq d<P ~ ^/Stt^ Jq d<P^ - Q" 

where for a system of finite total mass Qsup = 0 for a > 1, and Qsup = oo for a ^ 1 
(see the following discussion). Here (P is the (self-consistent) gravitational potential, Q = 
E + J^/2r2, where E and J are the specific (per unit mass) particle energy and angular 
momentum modulus, respectively, 

Pair) = + p(r) (2.2) 

is the so-called augmented density, and is the anisotropy radius (Osipkov 1979; 
Merritt 1985). Due to the dependence of / on global integrals of motion, it follows 
from the Jeans theorem that the equilibrium phase-space distribution function fiQ) 
satisfies the Vlasov equation (also known as the Collisionless Boltzmann Equation, see 
e.g. Binney & Tremaine 2008, Bertin 2014) Integration of fiQ) over the velocity space 
after multiplication by the appropriate velocity components shows that the anisotropy 
profile associated with eq. (2.1) is 


/3(r) = 1 - 


rrUr) 

2cr2(p) 



(2.3) 


where cr((r) and arir) are the tangential and radial velocity dispersion profiles of the 
system (see e.g. Binney & Tremaine 2008). In practice, Osipkov-Merritt models are 
everywhere isotropic for Ca —t oo and become more and more radially anisotropic for 
decreasing values of ra', at fixed they are isotropic for r <^ra and radially anisotropic 
for r ;g> r^. It is easy to show that the above formulae, derived in the literature for 
Newtonian gravity, are independent of the specific force law adopted. 

From eq. (2.1) it follows that the first step needed to recover / is the knowledge of 
the self consistent potential ‘P of the system. For a generic density p(x) the a—potential 
can be written as 

<l>(x) - <I>(xo) =-^ ^^^p(y)clV (2.4) 

for a 7 ^ 1, and as 

<I>(x) - <I>(xo) = G / In I^i^Zip(y)d3y (2.5) 

Jk? l|xo-y|| 

for a = 1. In the equations above G is a dimensional coupling constant and xq is 
a reference point to be chosen from convergence considerations. The value of ??(xo) 
is usually fixed by arguments of convenience as illustrated below. Note that, when 
convergence is assured, eq. (2.5) follows from (2.4) for a —>■ 1, provided that ||xo|| < oo. 

Elementary integration shows that for a spherical density distribution p(r) of finite 
mass, the potential can be written for 1 < a < 3 as 


^(r) 


27rG 

r 



(r -I- r')3 “ — |r — r' 
(a — 1)(3 — a) 


,/|3-a 


-r' pir')dr', 


( 2 . 6 ) 


having set to zero the value of the potential for r —>■ oo, and assumed ||xo|| —>■ oo. For 
a = 1 some care is needed in order to define d>ir), and it can be shown that for a system 
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Figure 1. Phase-space distribution functions for Osipkov-Merritt radially anisotropic Hernquist 
models with a = 1, 1.5, 2 and 2.5. The normalization constant is Q* = where r, is the 

half-mass radius. The black curves correspond to the isotropic case (C — l)i while the green ones 
to ^ = 1.5, the hducial value for the onset of ROI in Newtonian (a = 2) systems. The red curves 
show f{Q) for values of ^ near the consistency limit = f{rac) (see Sect. 2.2). 


of finite mass the potential can be written as 


<l>{r) 



/ /x2, ^ + 

(r + r y In- 

Tn 


(r 



2rr' 


r'p(r')dr', 


(2.7) 


where now Xq = 0, and is an arbitrary normalization length: different choices of 
just correspond to a different additive constant in the value of As expected eqs. (2.6) 
and (2.7) are more complicated than their Newtonian counterparts, because Newton’s 
theorems do not hold for a y 2. However, in general the potential at large radii for 
1 < a < 3 is dominated by the very simple asymptotic monopole term 


<?(r) 


GM 

(a — l)r““3 ’ 


( 2 . 8 ) 


while in the a = 1 case the monopole term for r —>■ oo diverges as 

<?(r) ~ GM In—. 

r-n 


(2.9) 


Note that, the expansions in eqs. (2.8) and (2.9) hold true also in the case of non-spherical 
density profiles, as can be seen from eqs. (2.4) and (2.5) for systems of finite total mass. 

In the present work we focus on systems described by the Hernquist (1990) density 
distribution 


p(r) 


M 

27rr^r(l -|- r/rcY' 


( 2 . 10 ) 


where M is the total mass, Tc is the core radius and r* = (1 -|- •\/2)i'c is the half-mass 
radius, i.e. the radius of the sphere enclosing half of the total mass. 

For the distribution (2.10) and 1 < a < 3 the potential at the centre can be evaluated 
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analytically and reads 


GM2B{3-a,a) 
rc~^ a — 1 


( 2 . 11 ) 


where B(x,j/) is the Complete Euler Beta function. Curiously, the dimensionless factor 
in equation above is a non monotonic function of a, with maximum -1 reached for a = 2, 
and diverging at —oo for a —)■ 1’'" and a —)■ 3“. For a = 1 the central potential is instead 
finite and reads 


<P(0) = GM + (2.12) 

its numerical value depends on the normalization length r„. For Hernquist models with 
a = 1 the natural choice is to adopt r„ = Tc- It is not surprising that the limit of eqs. 
(2.8) and (2.11) for a —)■ 1+ differs from eqs. (2.9) and (2.12), as two different choices 
of Xg have been adopted for 1 < a < 3 and a = 1 (see the discussion after eq. 2.5). In 
general, the Chebyshev (1853) theorem on the integration of differential binomials assures 
that the integral in eq. (2.6) can be evaluated in terms of elementary functions for the 
Hernquist model for all rational values of a. However, the formulae are cumbersome, 
so we computed eqs. (2.6) and (2.7) numerically; as a safety check we verified that the 
numerical results reproduce the asymptotic expansions at large radii (equations 2.8 and 
2.9) as well as the values of <?(0). 


2.2. Consistency 

In Newtonian gravity (a = 2) anisotropic Osipkov-Merritt systems of given p{r) are 
characterized by a critical value of the anisotropy radius, Vac, such that for < Vac 
the corresponding models are inconsistent (Ciotti & Pellegrini 1992, see also Carollo, de 
Zeeuw & van der Marel 1995; Ciotti 1996, 1999; Ciotti & Morganti 2010a,b; An, van Hese 
& Baes 2012). It is therefore important for our study to determine the minimum value 
Tacict) associated with a nowhere negative (i.e. consistent) f{Q). From eq. (2.1) it follows 
that the phase-space distribution function of an Osipkov-Merritt model (independently 
of the force law) can be written as 

fiQ) = MQ) + ^^, (2.13) 


where fi is the phase-space distribution function in the isotropic case (r^ —>■ oo). For 
the Hernquist density profile we found that fi^O for all the explored values of a i.e. 
the isotropic Hernquist model is consistent under r~°‘ forces, as already known for the 
Newtonian case. Being this assured, then the critical anisotropy radius is given by 


rlc = sup 


' fa{Qy 
. MQ)\' 


QgA, 


(2.14) 


where the sup in the r.h.s. is evaluated over the region A~ of Q values corresponding to 
/a < 0 (see Ciotti 2000). 

The numerically recovered phase-space distribution functions for a = 1, 1.5, 2 and 2.5 
and different values of ra (with the corresponding values of ^) are presented in Fig. 1. 
Of course, <?( 0 ) < Q < 0 for 1 < a < 3, and ^( 0 ) < Q < oo for a = 1. For ^ ~ be. 
Ta — Vac (red curves), f(Q) shows a strong non-monotonicity, due to the almost dominant 
negative contribution of fa/fa over the isotropic and positive fi. A further reduction of 
Ta would deepen even more the depression in the red curves, finally leading to a negative 
/, and to phase-space inconsistency. The relation of monotonicity of f{Q) and ROI is a 
very important point that at this stage we cannot address, but that will be central in 
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Figure 2. Minimum value r^c of the anisotropy radius for phase-space consistency (in units of 
the half-mass radius r,) of the Hernquist model with Osipkov-Merritt radial anisotropy and the 
corresponding value of the Fridman-Polyachenko-Shukhman index = C(i’ac), as functions of 
the force exponent a. 
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Figure 3. Fridman-Polyachenko-Shukhman parameter ^ as a function of the normalized 
anisotropy radius ra/r^, and of the force exponent a, for Osipkov-Merritt radially anisotropic 
Hernquist models with forces. 


our forthcoming study (for example by considering the possibility to extend the Antonov 
laws to r““ forces; see Binney & Tremaine 2008). 

In Fig. 2 we show for 1 < a < 3 the values of rac, and their associated = f{fac)- 
In the Newtonian case (a = 2) radT* — 0.08, in accordance with analytic results (Ciotti 
1996). From Fig. 2 it is also apparent that systems with lower values of a (i.e. forces 
with a longer range than Newtonian) can sustain smaller values of the anisotropy radius 
(associated with larger values of ^c) implying that they are able to store a larger fraction 
of the kinetic energy in low-angular momentum orbits maintaining a nowhere negative 
phase-space distribution function. 

Figure 3 (left panel) shows the parameter ^ as a function of the anisotropy radius 
for a = 1, 1.5, 2 and 2.5. As expected, ^ decreases monotonically for increasing r^ for all 
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values of a. A complementary illustration of this behaviour is given in the right panel, 
where the trend of ^ with a for a selection of values of is given. 

A more detailed representation of the effects of a and on the internal dynamics 
of the models can be obtained by plotting the velocity dispersion profiles. The radial 
velocity dispersion is given by the Jeans equation 


dpcr^ ^ 2j3pa‘^ d^ 
dr r ^ dr ’ 

whose solution in case of Osipkov-Merritt anisotropic systems is 


p{r)a^{r) 


rPir')ir'^+r^a)lpdr' 


(2.15) 


(2.16) 


It is easy to show that, for a density profile decreasing at large radii as r ^ with A > 3, 
and in presence of r““ forces. 


?(r) 


GM / r3-“ 

r2+r2 \a + X-3 


a + A — 1 


(2.17) 


In particular, for a = 1, in the limit r —>■ oo, flattens to the constant value GM/{X — 2) 
due to the logarithmic nature of the potential for a system of Hnite mass (see eq. 2.9). 
Specifically, for the Hernquist density profile given by eq. (2.10), A = 4, and thus cr^ — 
a/ GM/2 at large radii. In Fig. 4 we present the radial and tangential velocity dispersion 
profiles ar{r) and at{r) for r^/r* = 0.5, 1 and 2, obtained by numerically solving eq. 
(2.16). As a further test of the numerical recovery of f{Q) we verified that the profiles 
of ar and at extracted from body realizations obtained by sampling the numerically 
evaluated f{Q) reproduce the curves shown in Fig. 4. It is apparent how for increasing 
a, independently of the components of the velocity dispersion become more and more 
peaked at inner radii. As expected, for the MOND-like (a = 1) case, the radial profile of 
ar is always flat for all the explored values of in agreement with eq. (2.17) and with 
the numerical results of Di Cintio, Ciotti & Nipoti (2013). Once the profiles of a^ and af 
are known (see eq. 2.3), it is immediate to evaluate Kr and Kt by numerical integration. 
In particular, the relatively flat profiles of ar at large radii for small values of a are 
at the origin of the high values of ^ with respect to those obtained for larger values of 
a. The behaviour of ^ as a function of ra and a can be also easily illustrated with the 
aid of a very simple toy model. In practice, one can integrate eq. (2.17) for a constant 
density ’’atmosphere”, p(r) = po for r ^ R, under the gravitational field of a central point 
mass, —GM/r°‘. The radial velocity dispersion profile is of elementary evaluation, while 
Kr in general is given by hypergeometric (Euler Beta) functions. However, again from 
the Chebyshev theorem, Kr can be expressed in terms of elementary functions for all 
rational values of a. Kt is then obtained from the Virial theorem as Kt = —Kr + |IF|/2 
where W is given in the following eq. (3.2), which for the toy model here considered gives 
W = -47rGpoMR^-°‘/(4 - a). 


3. Stability of isotropic systems 

As a first sample of the numerical simulations that will be performed in a forthcoming 
work, we present here direct A—body simulations aimed at testing the stability of 
isotropic (^ = 1) Hernquist models for various values of a. The initial conditions are 
generated as in Nipoti, Londrillo & Ciotti (2002) with a Monte Carlo sampling of f{Q). 
The equations of motion for A = 25000 equal-mass particles are integrated with a second 
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ra/r*=0.5 ra/r»=1.0 ra/r»=2.0 



10-1 ;^oO ]^ol ;^o2 

r/r» 



Figure 4. Radial and tangential velocity dispersion profiles of the self-consistent radially 
anisotropic Hernquist model with forces, for ra/r* = 0.5, 1 and 2; r* is the half-mass 
radius of the distribution. The velocity dispersions are expressed in units of u* = r*/tt, and t, 
is given in eq. (3.1). The flat outer parts of the green curve (a = 1) in the top panels are in 
excellent agreement with the expected value from eq. (2.17), and A = 4. 


order symplectic scheme with fixed timestep At = t*/100, where the natural dynamical 
time-scale is defined as 


— 



(3.1) 


(see Di Cintio & Ciotti 2011 and Di Cintio, Ciotti & Nipoti 2013). In the code, the 
divergence of force and potential at vanishing inter-particle separation is cured with the 
introduction of a softening length e, so that in practice the potential due to each particle 
is smoothed as ^ oc l/(r^ -|- e 2 ^(a-i )/2 q, > ^ ln(\/r^ -|- e^) for a = 1. The 

optimal e (in units of r*) is fixed so that the softened force on a particle placed at ~ 5r* 
from the centre of mass of the system differs less than 0.01% from the unsoftened force. 
Therefore the softening length e is different for different a: e = 0.05, 0.07, 0.1 and 0.125 
for 0 = 1, 1.5, 2 and 2.5 respectively. 

The absence of the analogue of the Poisson equationf for r““ forces and a 2, 
points at direct TV— body simulations as the most natural approach for numerical studies. 
However, it is also important to note the well known property that potentials associated 
with the force can be expanded in terms of Gegenbauer polynomials for a 1, and 
of Chebyshev polynomials for a = 1 (Arfken & Weber 2005), thus opening the way to 
the construction of tree-codes (see e.g., Srinivasan et al. 2005). 

In the upper panel of Fig. 5 we show, for some representative value of a, the time 
evolution of the minimum-to-maximum semiaxis ratio c/a of the systems. The semiaxes 
are defined from the eigenvalues h ^ h h of the tensor lij = J2k calculated 

within the radius containing 0.85% of the total mass M (Meza & Zamorano 1997; Nipoti, 
Londrillo & Ciotti 2002). In the lower panel of the same figure we also show the time 
evolution of the corresponding virial ratio 2iV/|lT|, where K is the total kinetic energy 


f More specifically r “ forces do not obey a Poisson-like equation, belonging to the class of 
the so-called Riesz potentials, associated to fractional Laplace operators (Stein 1970). 
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0 5 10 15 20 25 30 35 40 


t/t. 

Figure 5. Evolution of the fiducial axial ratio c/a and the virial ratio 2K/\W\ for body 
simulations of isotropic Hernquist models with forces, and a = 1, 1.5, 2, 2.5. 


and 


W = 


N 

p(x)(x, V^>)d^x = ^ mi{xi,aji) 


(3.2) 


is the virial function. The second expression holds for a discrete system of particles with 
masses rrii, where 


- —GfTlj 



(3.3) 


is the acceleration on particle i due to particle j. We recall that for a 1 IT is related 
to the system total potential energy U by the identity 


W = {a- 1)17, 


(3.4) 


where 



1 


(3.5) 


and in the case of a discrete system <Pji is the potential on particle i due to particle j. For 
the a = 1 case (as for systems in deep-MOND regime) the virial function is independent 
of time, being: 


for a continuous system, and 



(3.6) 


W = - 


G 

~2 


n 

mi-nij 


(3.7) 


for a discrete system (see Di Cintio, Ciotti & Nipoti 2013). 

From Fig. 5 it is apparent that the spherical symmetry as well as the equilibrium of the 
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Figure 6. Final (t/t* = 40) density profiles of the same systems as in Fig. 5. Here p* = p(r*). 
The black solid line is the initial density profile given by eq. (2.10). The vertical dashed lines 
indicate the softening length e used in the simulations. 


initial conditions is preserved for all the considered values of a up to = 40, when the 
simulations end. The fluctuations are of the order of ~ 5% in c/a and ~ 1% in 2Kl\yV\. 
The stability of the isotropic systems is confirmed also by Fig. 6 where we compare the 
final density profile with the initial profile given by eq. (2.10). The initial density profile is 
preserved down to r/r» ~ 0.025 in the best case (a = 1.5) and to r/r* ~ 0.8 in the worst 
case {a = 2.5). In all cases the final density profiles are indistinguishable from the initial 
profiles for r > 1.4e. In the forthcoming study we will perform a quantitative monitoring 
of the deterioration of the central profile as a function of a, e and fV, by considering the 
time evolution of the Lagrangian radii containing a fixed fraction of the total mass. 


4. Discussion and conclusions 

In this preparatory exploration of the ROI in spherically symmetric collisionless 
equilibrium systems with forces we limited ourselves to a first important step, 
i.e. the construction of Osipkov-Merritt radially anisotropic self-consistent equilibrium 
models. We have recovered the potential generated by the Hernquist density distribution 
under forces and we have calculated the corresponding phase-space distribution 
functions for isotropic and anisotropic Osipkov-Merritt models in the range of values 
1 ^ a < 3. We numerically determined, as a function of a, the minimum value of the 
anisotropy radius for consistency rac and the associated critical value of the Fridman- 
Polyachenko-Shukhman index and we computed the radial and tangential velocity 
dispersion profiles for a few representative consistent systems. The main results of this 
paper can be summarized as follows: 

(i) For 1 ^ a < 3 the distribution functions of the isotropic Hernquist models 
are monotonically decreasing and everywhere positive. Their body realizations are 
numerically stable up to ~ 40f*. 

(ii) For anisotropic Hernquist models the minimum value of the anisotropy radius 
for consistency r^c increases with increasing a, and the corresponding fc decreases. 
Therefore, systems with lower a can be generated with a higher amount of radial 
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anisotropy. This appears to be in agreement with the results of Di Cintio, Ciotti & Nipoti 
(2013), who find increasingly anisotropic end products of dissipationelss collapses, when 
reducing a from 3 to 1. 

(iii) As general trend, for all values of a the Fridman-Polyachenko-Shukhman index ^ 
decreases for increasing ra, and increases with decreasing a at fixed ra¬ 
in the natural follow-up of this preparatory work we will determine, as a function of a, 
the critical ^ for the onset of ROI, in order to understand how much this process depends 
on the short or long range nature of the interaction law and we will study the structural 
and dynamical properties of the final systems. 
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